; docformat = 'rst'
;
; NAME:
;   MR_FIND_COUNTRY
;
; PURPOSE:
;   Find the country that resides in each cell in a specified grid
;
;+
; :Categories:
;   Model utilities
;
; :Returns:
;   2D array of country indices.  The names corresponding to each country can be found
;   by setting the Name keyword
;
; :Params:
;   LonIn : longitude of center for each grid cell in degrees
;   LatIn : latitude of center for each grid cell in degrees
;
; :Keywords:
;   Name : (output) Name of country corresponding to each returned index
;   Continent : (input) set this to TRUE to return continents instead of countries
;   Ocean : (input) set this to TRUE to allocate ocean gridcells to nearest country
;
; :Requirements::
;   Requires mr_country.gif, mr_continent.gif, mr_country_info.sav, mr_country_ocean.gif files
;   
;   Regular grids only.
;
; :Examples:
;   IDL> country_index=mr_find_country(lon, lat)
;
; :History:
;   Written by: Matt Rigby, MIT, May 29th 2011
;   With thanks to the GAMAP developers, whose map I've stolen!
;     http://acmg.seas.harvard.edu/gamap/
;
;-
function mr_find_country, lonin, latin, name=name, continent=continent, ocean=ocean, res=res

  ;Find location of this routine to find data files
  Help, Calls=callStack
  callingRoutine = (StrSplit(StrCompress(callStack[0])," ", /Extract))[0]
  thisRoutine = Routine_Info(callingRoutine, /Functions, /Source)
  thisDirectory=strmid(thisRoutine.path, 0, strpos(thisRoutine.path, '/', /reverse_search)+1)

  if keyword_set(res) then begin
    fname_res='_' + res
  endif else begin
    fname_res=''
  endelse

  if keyword_set(ocean) then ocean_country=1 else ocean_country=0 ;To avoid overwriting keyword with restored OCEAN variable
  if keyword_set(continent) then if continent eq 1 then continent_output=1 else continent_output=0 else continent_output=0 

  lon=lonIn
  Lat=latIn

  LonSize=n_elements(lon)
  LatSize=n_elements(lat)
  dlon=lon[1] - lon[0]  ;ONLY REGULAR GRIDS
  dlat=lat[1] - lat[0]
  
  ;Read data files
  restore, strcompress(thisDirectory + 'mr_country_info' + fname_res + '.sav')
  read_gif, strcompress(thisDirectory + 'mr_country' + fname_res + '.gif'), country
  country=fix(country)
  read_gif, strcompress(thisDirectory + 'mr_country_ocean' + fname_res + '.gif'), ocean
  ocean=fix(ocean)

  ;make country longitude match input longitude
  if total(lon lt 0.) then begin
    lon_temp=lonc
    lat_temp=latc
    mr_grid_orient, lon_temp, lat_temp, country, /neglon
    mr_grid_orient, lonc, latc, ocean, /neglon
  endif


  hlon=histogram(lonc, min=(lon[0]-dlon/2.), binsize=dlon, nBins=n_elements(lon), reverse_indices=lonIndex)
  hlat=histogram(latc, min=(lat[0]-dlat/2.), binsize=dlat, nBins=n_elements(lat), reverse_indices=latIndex)

  country_out=fltarr(LonSize, LatSize)
  ocean_out=fltarr(LonSize, LatSize)

  for LonI=0, LonSize-1 do begin
  for LatI=0, LatSize-1 do begin
    if (lonIndex[lonI] ne LonIndex[LonI+1]) and (latIndex[latI] ne LatIndex[LatI+1]) then begin
      LonII=LonIndex[lonIndex[lonI]:LonIndex[LonI+1]-1]
      LatII=LatIndex[latIndex[latI]:LatIndex[LatI+1]-1]
      country_out[LonI, LatI]=mr_mode(country[lonII[0]:LonII[-1],latII[0]:LatII[-1]])
      ocean_out[LonI, LatI]=mr_mode(ocean[lonII[0]:LonII[-1],latII[0]:LatII[-1]])
    endif
  endfor
  endfor

  wh=where(country_out eq 255, count)
  if ocean_country then begin
    if count gt 0 then begin
      country_out[wh]=ocean_out[wh]
    endif
  endif
  index=country_out[uniq(country_out, sort(country_out))]
  index=index[where(index ne 255)]

  name=countryname[index]
  continent=continentname[continentlookup[index]]
  country_out[where(country_out eq 255)]=!values.f_nan

  ;re-order country number
  for i=0, n_elements(index)-1 do country_out[where(country_out eq index[i])]=i
  index=indgen(n_elements(index))


  ;OUTPUT CONTINENT NUMBER in COUNTRY_OUT and CONTINENT NAME IN 'NAME'
  if continent_output then begin
    continent_out=replicate(255., lonsize, latsize)
    for ci=0, n_elements(index)-1 do begin
      wh=where(country_out eq index[ci])
      if wh[0] ne -1 then begin
        continent_out[wh]=where(continentname eq continent[ci])
      endif else print, 'WHAT?!', ci
    endfor

    index=continent_out[uniq(continent_out, sort(continent_out))]
    index=index[where(index ne 255)]
    continent_out[where(continent_out eq 255)]=!values.f_nan
    name=continentname[index]

    ;re-order country number
    for i=0, n_elements(index)-1 do continent_out[where(continent_out eq index[i])]=i
    index=indgen(n_elements(index))

    country_out=continent_out
  endif


  return, country_out

end